function [ results_fnames ] = osl_run_first_level_continuous_state_vector_seed( oat )

% [ results_fnames ] = osl_run_first_level_continuous( oat )
%
% takes in an OAT, which needs to be setup by calling oat=osl_setup_oat(S), struct
% and runs first level continuous time GLM or seed based connectivity
%
% This function should normally be called using osl_run_oat(oat);
%
% MWW 2012

global OSLDIR;

do_glm_statewise=oat.first_level.hmm_do_glm_statewise;

first_level=oat.first_level;
source_recon_name='source_recon';

if strcmp(oat.source_recon.modalities{1},'EEG') 
    modality_meeg='EEG';
else
    modality_meeg='MEG';
end

% check if user has asked not to do GLM
if(~isfield(first_level,'doGLM')),
    first_level.doGLM = 1;
end

if do_glm_statewise,
    error('Doing do_glm_statewise is not implemented with seed-based connectivity');
end;

if(first_level.time_average)
    error('Time averaging over whole window is not appropriate on continuous data');   
end;

if(~isfield(first_level,'connectivity_seed_mni_coord')),
    error('connectivity_seed_mni_coord needs specifying');
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set the downsample factors 
% tf_downsample_factor is used in source space after recon weights are applied 
% time_downsample_factor is used after TF transform (as part of osl_tf_transform call) in sensor space on complex values
if ~isempty(first_level.post_tf_downsample_factor)
    post_tf_ds_factor = 1./first_level.post_tf_downsample_factor;
else
    post_tf_ds_factor = 1;
end

if ~isempty(first_level.post_movingaverage_downsample_factor)
    post_movingaverage_ds_factor = 1./first_level.post_movingaverage_downsample_factor;
else
    post_movingaverage_ds_factor = 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ntrials=1;

triallist=1;
   
clear stats;

for subi_todo=1:length(first_level.sessions_to_do),   
        
    sub=first_level.sessions_to_do(subi_todo);
    
    %sub/length(first_level.submatfiles),
    
    %% load in the beamformer result for this session
    source_recon_results_fname = oat.source_recon.results_fnames{sub};
    source_recon_results=osl_load_oat_results(oat,source_recon_results_fname);
    
    % get D and update path in case OAT dir has been moved
    D = osl_get_oat_sensordata( source_recon_results );
    
    %% setup things that are common to all sessions
    if(subi_todo==1),
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% set up mask indices from mask/coords
        if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
            
            error('Working in sensor space - not implemented for seed based connectivity');

            %%%%%%%%%%%%%%%
            %% setup channels
            chanindmeg = strmatch(modality_meeg, D.chantype);
            first_level_results.mask_indices_in_source_recon=1:length(chanindmeg);

        else,

            if(isfield(first_level,'mni_coords'))    
                                
                current_level_mni_coord = first_level.mni_coords;
                
                % load in lower level coords
                lower_level_mni_coord=source_recon_results.mni_coord;               
                %lower_level_mask_fname  = [oat.source_recon.dirname '/' source_recon_name '_mask'];
                %[ lower_level_mni_coord, ~ ] = osl_mnimask2mnicoords(lower_level_mask_fname);

                [~, iA, mask_indices]=intersect(current_level_mni_coord, lower_level_mni_coord,'rows');

                % sort indices to be in the order of the lower level mni_coords:
                [~, gg]=sort(iA);
                indices_in_lower_level=mask_indices(gg);

                first_level_results.mask_indices_in_source_recon = indices_in_lower_level;
                first_level_results.mni_coord = lower_level_mni_coord(indices_in_lower_level,:);
                
                % the mni coordinates from the specified list
                str=['Desired MNI coordinates ' ];
                for vox=1:size(current_level_mni_coord,1),
                    str=[str ', [' num2str(current_level_mni_coord(vox,:)) ']'];
                end;
                disp(str);
                
                % the mni coordinates from the lower level mask
                str=['Nearest dipoles from lower level at MNI coordinates ' ];
                for vox=1:length(indices_in_lower_level),
                    str=[str ', [' num2str(lower_level_mni_coord(indices_in_lower_level(vox),:)) ']'];
                end;
                disp(str);
            
            else
                % setup std space brain
                first_level_results.gridstep=source_recon_results.gridstep;

                S=[];
                S.lower_level_mask_fname=[oat.source_recon.dirname '/' source_recon_name '_mask'];
                S.current_level_mask_fname=[oat.source_recon.dirname '/' oat.first_level.name '_mask'];
                S.current_level=first_level;
                S.lower_level_mni_coord=source_recon_results.mni_coord;
                S.lower_level_gridstep=source_recon_results.gridstep;

                [first_level_results.mask_indices_in_source_recon, first_level_results.mni_coord]=setup_mask_indices(S);
                lower_level_mni_coord=source_recon_results.mni_coord;
                clear S;
            end;
        end;

        %%%%%%%%%%%%%%%%
        %% setup seed for seed-based connectivity analysis       
        if(isempty(first_level.connectivity_seed_mni_coord)),
           error('need to specify seed');
        end;
 
        mni_coords=lower_level_mni_coord;
        
        disp(['Using for seed: MNI coordinate ' num2str(first_level.connectivity_seed_mni_coord)]);

        % find the nearest index of the beamformed voxels to the specified mni_coord
        dists=(sqrt(sum((mni_coords-repmat(first_level.connectivity_seed_mni_coord,size(mni_coords,1),1)).^2,2)));

        [dist,first_level_results.seed_index]=min(dists);

        disp(['Nearest dipole for seed at MNI coordinate ' num2str(mni_coords(first_level_results.seed_index,:))]);                                              

        first_level.contrast=[];
        first_level.contrast{1}=[1]';

        %%%%%%%%%%%%%%%%
        %% results containers
        Nvoxels_out=length(first_level_results.mask_indices_in_source_recon);

        first_level_results.pseudo_zstat_var=zeros([Nvoxels_out,1]);
        
    end; % end of if(sub==1)

    disp(['Reconstruct time courses and computing stats for dataset ' source_recon_results.fname]);    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% setup HMM 
    if strcmp(source_recon_results.recon_method,'none'),
        NKrecon=1;
    else
        NKrecon=numel(source_recon_results.BF.inverse.W.MEG);                          
    end;

    if strcmp(source_recon_results.recon_method,'none'),
        NKrecon=1;
        NKglm=NKrecon;
    else
        if strcmp(oat.source_recon.modalities{1},'EEG')
            NKrecon=numel(source_recon_results.BF.inverse.W.EEG);
        else
            NKrecon=numel(source_recon_results.BF.inverse.W.MEG);
        end            

        NKglm=1;            
    end;
        

    %%%%%%%%%%%%%%%
    %% setup channels
    chanindmeg = strmatch(modality_meeg, D.chantype);
    chanind = setdiff(chanindmeg, D.badchannels);
    if isempty(chanind)
        error(['No good ' modality_meeg ' channels were found.']);
    end
    if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
        % for sensor space - also include bad sensors sensors 
        first_level_results.mask_indices_in_source_recon=1:length(chanindmeg);
        chanind = chanindmeg;
    end;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% set up potential frequency bins
    fsample=D.fsample;

    if(~isempty(first_level.tf_hilbert_freq_ranges)),
        disp('oat.first_level.tf_hilbert_freq_ranges is set, freq_range, tf_freq_res and tf_num_freqs will be overridden');

        first_level.tf_num_freqs=size(first_level.tf_hilbert_freq_ranges,1);
        freq_range=[]; 
    else
        freq_range = first_level.tf_freq_range; % range used for stat testing
        if isempty(freq_range)
            % if not set as a first_level option, default to the range used
            % for source_recon
            freq_range=source_recon_results.source_recon.freq_range;

            if isempty(freq_range)
                freq_range=[1 fsample/3];
                disp(['WARNING: No frequency range has been specified. Set oat.source_recon.freq_range or oat.first_level.tf_freq_range. Using ' num2str(freq_range) ' Hz.']);                
            end % isempty(freq_range)
        end % if isempty(freq_range)

    end;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% set up time windows
    % Time windows:
    % - source_recon_results.source_recon.time_range is all of the desired
    % time range for the source recon: this is used as the input time range into tf
    % transform
    % - first_level.time_range is all of the desired
    % time range for the first level.
    % - source_recon_results.samples2use and source_recon_results.woi 
    % should be the same, and are the time windows corresponding to only 
    % the good epochs within source_recon_results.source_recon.time_range.
    
    % For feeding into tf, we want to include even the bad epochs, i.e. use 
    % source_recon_results.source_recon.time_range (best compromise wrt
    % edge effects)
   
    if(isempty(source_recon_results.source_recon.time_range)),
        source_recon_results.source_recon.time_range=[D.time(1) D.time(end)];
    end;
    
    source_time_range=source_recon_results.source_recon.time_range;
   
    tims.D.times=D.time;
    tims.D.to_tf_input=( D.time >= source_time_range(1) ) & ( D.time <= source_time_range(2) ); % time indices in D.time for use by tf
    tims.tf_input.times=tims.D.times(tims.D.to_tf_input);
    
    disp(['Sample rate of sensor data is ' num2str(fsample) ' Hz']);
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% setup post tf time and frequency windows/indices for using later
    
    if isempty(first_level.time_range)
        first_level.time_range=source_recon_results.source_recon.time_range;
    end
        
    if(size(first_level.time_range,1)>1)
        error('first_level.time_range can only be a single time window');
    end;

    % Setup tf output time dimension.   
    % out.tf_times is time indices output from tf transform. These should
    % now be adjusted to exclude bad epochs and only work within the 
    % specified first_level.time_range.
    % So need to combine first_level.time_range and
    % source_recon_results.samples2use (source_recon_results.woi) to get
    % time indices of interest, and find intersect with out.tf_times
        
    % for testing: source_recon_results.woi=[250 350; 370 378; 400 800; 990 1000];
    
    % first combine source_recon_results.woi and first_level.time_range
    first_level_woi_tmp=source_recon_results.woi;
    first_level_woi=[];
    j=1;
    for i=1:size(first_level_woi_tmp,1)
        first_level_woi_tmp(i,1)=max(first_level_woi_tmp(i,1),first_level.time_range(1));
        first_level_woi_tmp(i,2)=min(first_level_woi_tmp(i,2),first_level.time_range(2));
        if(first_level_woi_tmp(i,2)>first_level_woi_tmp(i,1))
            first_level_woi(j,:)=first_level_woi_tmp(i,:);
            j=j+1;
        end;
    end;    
    
    % apply this to D.times 
    tims.tf_input.to_tf_input_clean = false(1,numel(tims.tf_input.times));
    for i=1:size(first_level_woi,1)
         tims.tf_input.to_tf_input_clean(( tims.tf_input.times >= first_level_woi(i,1) ) & ( tims.tf_input.times <= first_level_woi(i,2) )) = true;
    end
    tims.tf_input_clean.times=tims.tf_input.times(tims.tf_input.to_tf_input_clean);

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% set up tf basis

    tf_settings = [];
    tf_settings.tf_method               = first_level.tf_method;
    tf_settings.tf_logtransform         = first_level.tf_logtransform;
    tf_settings.tf_hilbert_freq_res     = first_level.tf_hilbert_freq_res;
    tf_settings.tf_hilbert_freq_ranges  = first_level.tf_hilbert_freq_ranges;
    tf_settings.tf_freq_range           = freq_range; % set above
    tf_settings.tf_num_freqs            = first_level.tf_num_freqs;
    tf_settings.tf_multitaper_ncycles   = first_level.tf_multitaper_ncycles;
    tf_settings.tf_hilbert_freq_ranges  = first_level.tf_hilbert_freq_ranges;
    tf_settings.raw_times               = tims.tf_input_clean.times;
    tf_settings.ds_factor               = post_tf_ds_factor; 
    tf_settings.tf_morlet_factor        = first_level.tf_morlet_factor;
    tf_settings.tf_multitaper_taper                = first_level.tf_multitaper_taper;
    tf_settings.tf_multitaper_twin                 = first_level.tf_multitaper_twin;
    tf_settings.tf_multitaper_freqsmooth           = first_level.tf_multitaper_freqsmooth;
    tf_settings.tf_calc_amplitude=0; % we will calc the amplitude after applying recon weights        
    tf_settings.tf_hilbert_do_bandpass_for_single_freq = first_level.tf_hilbert_do_bandpass_for_single_freq;

    out = osl_tf_transform( tf_settings ); % passed without a data matrix, osl_tf_transform just gets container sizes

    if isfield(out, 'tf_morlet_basis')
        first_level.tf_morlet_basis = out.tf_morlet_basis;
    end
    
    tims.tf_output_clean.times=out.tf_times;
       
    fsample_posttf=1/mode(diff(out.tf_times));

    disp(['Sample rate after TF transform is ' num2str(fsample_posttf) ' Hz']);

    timtims=tims.D.times;mat2str([min(timtims),max(timtims)])
    timtims=tims.tf_input.times;mat2str([min(timtims),max(timtims)])
    timtims=tims.tf_input_clean.times;mat2str([min(timtims),max(timtims)])
    timtims=tims.tf_output_clean.times;mat2str([min(timtims),max(timtims)])
       
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% averaging
    
    first_level_results.times=1;
    if first_level.time_average
        error('Time averaging over whole window is not appropriate on continuous data');   
    end

    % setup result freq dimension
    first_level_results.frequencies = out.tf_freqs;
    if first_level.freq_average,
        first_level_results.frequencies=mean(first_level_results.frequencies);
    end

    clear out;
        
    ntpts=1;
    nfreqs=length(first_level_results.frequencies); 

    %%%%%%%%%%%%%%%
    %% setup class windows
    classchanind=find(strcmp(D.chanlabels,'Class'));
    if(isempty(classchanind)),
        error(['No ''Class'' chanlabel in: ' D.fname]);
    end;
    
    % class_samples_inds_recon is used for using the appropriate weights in recon 
    for si=1:NKrecon,        
        class_samples_inds_recon_tf_input{si} = (D(classchanind, find(tims.D.to_tf_input), triallist)==si);  % has had bad epochs removed
    end;
 
    % class_samples_inds_glm is used for doing GLM. Can not do statewise
    % GLM - it is not implemented
    class_samples_inds_glm_tf_input{1}=ones(1,length(find(tims.D.to_tf_input)),length(triallist));
    class_samples_inds_glm_tf_input_clean{1}=class_samples_inds_glm_tf_input{si}(:,find(tims.tf_input.to_tf_input_clean),:,:);

    for si=1:NKrecon,
        % create clean (bad epochs removed) version for use later
        class_samples_inds_recon_tf_input_clean{si}=class_samples_inds_recon_tf_input{si}(:,find(tims.tf_input.to_tf_input_clean),:,:);

        state_tinds=class_samples_inds_recon_tf_input{si};                                 

        total_time=size(state_tinds,2)/fsample;
        state_time=total_time*sum(sum(squeeze(state_tinds)))/prod(size(squeeze(state_tinds)));
        txt=['State ' num2str(si) ' is active for ' num2str(state_time) 'secs'];

        if(state_time<10),
            warning([txt '. It will be excluded from the GLM analysis.']);
        else
            disp(txt);
        end;
    end;
    clear state_tinds;
    
    %%%%%%%%%%%%%%%
    %% load in sensor data
    first_level_results.D_sensor_data=D;
    first_level_results.chanind=chanind;
    
    sensor_data_tf_input=D(chanind, find(tims.D.to_tf_input), triallist); % note: still has bad epochs in it
    
    sensor_data_tf_input_clean=sensor_data_tf_input(:,find(tims.tf_input.to_tf_input_clean),:); 
    
    clear sensor_data_tf_input;
    
    %%%%%%%%%%%%%%%%%%%%%%
    %% recon the seed for doing orthogonalisation  
    if(isempty(first_level.connectivity_seed_mni_coord)),
        error('No seed specified');
    end;
        
    % recon seed time course
    S2=[];
    S2.source_recon_results=source_recon_results;
    S2.chanind=chanind;
    S2.class_samples_inds=class_samples_inds_recon_tf_input_clean;
    S2.voxind=first_level_results.seed_index;
    S2.reduce_rank=2;
    S2.num_freqs=1;
        
    Ntrials=size(class_samples_inds_recon_tf_input_clean{1},3);
    for kk=1:NKrecon,  
        for tri=1:Ntrials,
            S2.timeinds{kk,tri}=find(class_samples_inds_recon_tf_input_clean{kk}(1,:, tri)); % time indices for class kk         
            S2.sensor_data_sub{kk,tri}=sensor_data_tf_input_clean(:,S2.timeinds{kk,tri},tri);
        end;
    end;
    [ dat_seed_tf_input_clean wnorms_seed wnorm_seed wnorms_nai_seed wnorm_nai_seed, weights_seed ] = get_voxel_recon_timecourse_vector( S2 );
    
    clear S2;
    
    %dat_seed=dat_seed';
    %wnorms_nai_seed=wnorms_nai_seed';

    dat_seed_tf_input_clean=dat_seed_tf_input_clean./sqrt(nanmean(squash(wnorms_nai_seed))); % trials x tpts x dipole_rank
        
    %conn.dat_seed=permute(dat_seed_tf_input_clean(:,find(tims.tf_input.to_tf_input_clean),:,:),[2 3 1 4]); % clean: has bad epochs removed
    conn.dat_seed=permute(dat_seed_tf_input_clean,[2 3 1 4]); % clean: has bad epochs removed
    
    % tf seed in source space:
    tf_settings = [];
    tf_settings.tf_method               = first_level.tf_method;
    tf_settings.tf_logtransform         = first_level.tf_logtransform;
    tf_settings.tf_hilbert_freq_res     = first_level.tf_hilbert_freq_res;
    tf_settings.tf_hilbert_freq_ranges  = first_level.tf_hilbert_freq_ranges;
    tf_settings.tf_freq_range           = freq_range; % set above
    tf_settings.tf_num_freqs            = first_level.tf_num_freqs;
    tf_settings.raw_times               = tims.tf_input_clean.times;
    tf_settings.ds_factor               = post_tf_ds_factor; 
    tf_settings.tf_morlet_factor        = first_level.tf_morlet_factor;
    tf_settings.tf_multitaper_taper                = first_level.tf_multitaper_taper;
    tf_settings.tf_multitaper_twin                 = first_level.tf_multitaper_twin;
    tf_settings.tf_multitaper_freqsmooth           = first_level.tf_multitaper_freqsmooth;
    tf_settings.tf_calc_amplitude=0;
    tf_settings.tf_hilbert_do_bandpass_for_single_freq = first_level.tf_hilbert_do_bandpass_for_single_freq;
    
    if isfield(first_level,'tf_morlet_basis')
        tf_settings.tf_morlet_basis = first_level.tf_morlet_basis;
    end

    dims=size(dat_seed_tf_input_clean);
    %dims(2)=length(tims.tf_output_clean.times);
    %dat_seed_tf_output_clean=zeros(dims);
    clear dat_seed_tf_output_clean;
    for tri=1:size(dat_seed_tf_input_clean,3),
        for rnk=1:size(dat_seed_tf_input_clean,4),
            out = osl_tf_transform( tf_settings , permute(dat_seed_tf_input_clean(1,:,tri,rnk), [1,2,3,4]) ); % do the transformation
            dat_seed_tf_output_clean(tri,:,:,rnk)=permute(out.dattf,[2 1 3 4]);
        end;
    end;
              
    %conn.dat_seed_tf=abs(permute(dat_seed_tf_output(1,find(tims.tf_output.to_tf_output_clean),:,:),[2 3 1 4])); % cleaned - bad epochs removed
  
    conn.dat_seed_tf=abs(permute(dat_seed_tf_output_clean(1,:,:,:),[2 3 1 4])); % cleaned - bad epochs removed
  
    %%%%%%%%%%%%%%%
    %% Orthogonalise sensor space data wrt seed:
    sensor_data_tf_input_clean_zerolagcorr=zeros(size(sensor_data_tf_input_clean));
    if(first_level.connectivity_seed_regress_zerolag),

        do_seed_regress_zerolag_statewise=0;
        
        % orthoganalise the target voxel WRT the seed voxel (doing
        % it this way round means that the TF transform and
        % windowing does not need to be redone for the seed voxel)
        
        if(NKrecon>1 && do_seed_regress_zerolag_statewise),  
            
            for si=1:NKrecon,
                
                %tinds=(class_samples_inds_recon_tf_input{si}==1);
                tinds_clean=(class_samples_inds_recon_tf_input_clean{si}==1);
                
                for tri=1:size(sensor_data_tf_input_clean,3),
                
                    f=1; % pre-tf-transform

                    %dat2=permute(dat_seed_tf_input(tri,tinds,f,:),[2,4,3,1]); % ntpts x ndipoles
                    dat2_clean=permute(dat_seed_tf_input(1,tinds_clean,f,:),[2,4,3,1]); % ntpts x ndipoles

                    % iterate over sensors
                    for indind=1:size(sensor_data_tf_input_clean,1),
                        %dat1=permute(sensor_data_tf_input(indind,tinds,tri),[2,1,3]); % ntpts x 1
                        dat1_clean=permute(sensor_data_tf_input_clean(indind,tinds_clean,tri),[2,1,3]);

                        if(indind~=first_level_results.seed_index)

                            % iterate over dipole directions
                            for dd=1:size(dat2_clean,2),
                                % PE=dat2\dat1 carried out on clean data
                                % This then applied to unclean data to get the
                                % output
                                dat1_clean=dat1_clean-dat2_clean(:,dd)*(dat2_clean(:,dd)\dat1_clean);
                            end;

                        end;


                        sensor_data_tf_input_clean_zerolagcorr(indind,tinds_clean,tri)=dat1_clean;
                    end;
                
                end;
            end;
            
        else,
                 
            for tri=1:size(sensor_data_tf_input_clean,3),

                f=1; % pre-tf-transform

                % calc orthogonalisation regression params from cleaned data and apply
                % to non-clean
                %dat2=permute(dat_seed_tf_input(tri,:,f,:),[2,4,3,1]); % ntpts x ndipoles
                dat2_clean=permute(dat_seed_tf_input_clean(1,:,f,:),[2,4,3,1]); % ntpts x ndipoles

                % iterate over sensors
                for indind=1:size(sensor_data_tf_input_clean,1),
                    %dat1=permute(sensor_data_tf_input(indind,:,tri),[2,1,3]); % ntpts x 1
                    dat1_clean=permute(sensor_data_tf_input_clean(indind,:,tri),[2,1,3]);

                    if(indind~=first_level_results.seed_index)

                        % iterate over dipole directions
                        for dd=1:size(dat2_clean,2),
                            % PE=dat2\dat1 carried out on clean data
                            % This then applied to unclean data to get the
                            % output
                            dat1_clean=dat1_clean-dat2_clean(:,dd)*(dat2_clean(:,dd)\dat1_clean);
                        end;

                    end;

                    sensor_data_tf_input_clean_zerolagcorr(indind,:,tri)=dat1_clean;
                end;

            end;    
        end;
    end;        
    
    %%%%%%%%%%%%%%%
    %% Do T-F transform in sensor space
    %% Transform the data - either downsampling in time domain, or a TF transform
    
    disp('Doing T-F transform in sensor space...');

    tf_settings = [];
    tf_settings.tf_method               = first_level.tf_method;
    tf_settings.tf_logtransform         = first_level.tf_logtransform;
    tf_settings.tf_hilbert_freq_res     = first_level.tf_hilbert_freq_res;
    tf_settings.tf_hilbert_freq_ranges  = first_level.tf_hilbert_freq_ranges;
    tf_settings.tf_freq_range           = freq_range; % set above
    tf_settings.tf_num_freqs            = first_level.tf_num_freqs;
    tf_settings.raw_times               = tims.tf_input_clean.times;
    tf_settings.ds_factor               = post_tf_ds_factor; 
    tf_settings.tf_morlet_factor        = first_level.tf_morlet_factor;
    tf_settings.tf_multitaper_taper                = first_level.tf_multitaper_taper;
    tf_settings.tf_multitaper_twin                 = first_level.tf_multitaper_twin;
    tf_settings.tf_multitaper_freqsmooth           = first_level.tf_multitaper_freqsmooth;
    tf_settings.tf_calc_amplitude=0;
    tf_settings.tf_hilbert_do_bandpass_for_single_freq = first_level.tf_hilbert_do_bandpass_for_single_freq;
    
    if isfield(first_level,'tf_morlet_basis')
        tf_settings.tf_morlet_basis = first_level.tf_morlet_basis;
    end

    dims=size(sensor_data_tf_input_clean_zerolagcorr);
    sensor_data_tf_output_clean_zerolagcorr=zeros(dims(1),length(tims.tf_output_clean.times),Ntrials,first_level.tf_num_freqs);
    clear sensor_data_tf_output_clean_zerolagcorr;
    ft_progress('init', 'etf');
    for indind=1:size(sensor_data_tf_input_clean_zerolagcorr,1),            
        ft_progress(indind/size(sensor_data_tf_input_clean_zerolagcorr,1));

        for tri=1:size(sensor_data_tf_input_clean_zerolagcorr,3),
            out = osl_tf_transform( tf_settings , sensor_data_tf_input_clean_zerolagcorr(indind,:,tri) ); % do the transformation
            sensor_data_tf_output_clean_zerolagcorr(indind,:,tri,:,:)=permute(out.dattf,[2 3 1]);
        end;
    end;   
  
    ft_progress('close');    
    
    %%%%%%%%%%%%%%%%
    %% Temporally downsample state time courses (data will be done later after taking abs (for t-f data) inside the voxel loop)     
    if post_tf_ds_factor~=1 || post_movingaverage_ds_factor~=1,
        if(NKglm>1)
            error('Cannot have multiple classes and do downsampling');
        end;
   
        if(NKrecon>1)
            error('Cannot have multiple classes and do downsampling');
        end;
    end;
    class_samples_inds_glm_tf_output_clean{si}=ones(1,length(tims.tf_output_clean.times),Ntrials);
    class_samples_inds_recon_tf_output_clean{si}=ones(1,length(tims.tf_output_clean.times),Ntrials);
    
    %%%%%%%%%%%%%%%%
    %% do everything from now on one voxel at a time (for the sake of
    %% RAM) in the first level mask space                
    disp('Looping over voxels...');

    ft_progress('init', 'etf');

    % debug:
    %[indfind, mni_coord]=nearest_vec(source_recon_results.mni_coord,[-32 -26 50]);
           
    S2=[]; % S2 will be time domain, S3 with be TF space
    S2.source_recon_results=source_recon_results;
    S2.chanind=chanind;
    S2.class_samples_inds=class_samples_inds_recon_tf_input_clean;
    S2.reduce_rank=2;
    S2.num_freqs=first_level.tf_num_freqs;
        
    S3=S2;
    
    Ntrials=size(class_samples_inds_recon_tf_input_clean{1},3);
    for kk=1:NKrecon,  
        for tri=1:Ntrials,        
            S2.timeinds{kk,tri}=find(class_samples_inds_recon_tf_input_clean{kk}(1,:, tri)); % time indices for class kk
            S3.timeinds{kk,tri}=find(class_samples_inds_recon_tf_output_clean{kk}(1,:, tri)); % time indices for class kk
 
            S2.sensor_data_sub{kk,tri}=reshape(sensor_data_tf_input_clean_zerolagcorr(:,S2.timeinds{kk,tri},tri,:),[size(sensor_data_tf_input_clean_zerolagcorr,1),length(S2.timeinds{kk,tri})*size(sensor_data_tf_input_clean_zerolagcorr,4)]);
            S3.sensor_data_sub{kk,tri}=reshape(sensor_data_tf_output_clean_zerolagcorr(:,S3.timeinds{kk,tri},tri,:),[size(sensor_data_tf_output_clean_zerolagcorr,1),length(S3.timeinds{kk,tri})*size(sensor_data_tf_output_clean_zerolagcorr,4)]);

        end;
    end;
    
    need_non_tf_data=0;
    need_tf_data=1;
    
    if any(strcmp(oat.first_level.connectivity_fc_methods,'he')) || any(strcmp(oat.first_level.connectivity_fc_methods,'plv')) || any(strcmp(oat.first_level.connectivity_fc_methods,'synchro_hmm')) || any(strcmp(oat.first_level.connectivity_fc_methods,'hmm')),
    for indind=1:Nvoxels_out, % indexes brain space      
        
        ft_progress(indind/Nvoxels_out);                                   
        S2.voxind=first_level_results.mask_indices_in_source_recon(indind);

        if need_non_tf_data,
            [ dat wnorms wnorm wnorms_nai wnorm_nai weights] = get_voxel_recon_timecourse_vector( S2 );
        
            dat=dat./sqrt(nanmean(squash(wnorms_nai)));
        
            conn.dat=permute(dat(1,:,:,:),[2 3 1 4]);    
        end;
        
        S3.voxind=first_level_results.mask_indices_in_source_recon(indind);
        
        if need_tf_data,
            [ dat_tf wnorms wnorm wnorms_nai wnorm_nai weights] = get_voxel_recon_timecourse_vector( S3 );
            
            dat_tf=dat_tf./sqrt(mean(squash(wnorms_nai))); 

            conn.dat_tf=permute(dat_tf(1,:,:,:),[2 3 1 4]);
        end;        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
        %% now do the stats

        %%%%%%
        % on vec:

        tres=1/fsample_posttf;

        bandcentre=mean(first_level.tf_freq_range);

        if any(strcmp(oat.first_level.connectivity_fc_methods,'synchro_hmm')),
            %%
            bandcentre=mean(first_level.tf_freq_range);
            S=[];

            S.tres=tres;
            S.data_high=(conn.dat_seed);
            S.data_low=(conn.dat);
            S.NK=4;
            S.force_zero_means=1;
            S.num_hmm_starts=1;
            S.do_plots=0;
            S.hmm_pca_dim=10;                
            span=1/bandcentre;
            S.M=round(span/tres);
            S.deltasecs=span/S.M;

            %S.deltasecs=1/(S.M*bandcentre);

            res  = synchro_hmm( S );

            kl=zeros(S.NK);

            hmm=res.hmm;
            Apca=res.Apca;
            S.NK=S.NK;
            for ii=1:S.NK,     
                for jj=1:S.NK,
                    if(ii~=jj),
                        ci=hmm.state(ii).Cov;
                        cj=hmm.state(jj).Cov;
                        kl(ii,jj)=trace(pinv(cj)*ci)-log(det(ci)/det(cj))-size(cj,1);
                    end;
                end;
            end; 

            first_level_results.pw_oasm_embedded(indind)=sum(squash(kl))/(S.NK^2-S.NK);
            %%
        end;

        if any(strcmp(oat.first_level.connectivity_fc_methods,'hmm')),
            %%
            bandcentre=mean(first_level.tf_freq_range);
            S=[];

            S.tres=tres;
            S.data_high=(conn.dat_seed);
            S.data_low=(conn.dat);
            S.NK=4;
            S.force_zero_means=1;
            S.num_hmm_starts=1;
            S.do_plots=0;
            S.hmm_pca_dim=4;                
            S.deltasecs=tres;
            S.M=1;
            %S.deltasecs=1/(S.M*bandcentre);

            res  = synchro_hmm( S );

            kl=zeros(NK);

            hmm=res.hmm;
            Apca=res.Apca;
            NK=S.NK;
            for ii=1:NK,     
                for jj=1:NK,
                    if(ii~=jj),
                        ci=hmm.state(ii).Cov;
                        cj=hmm.state(jj).Cov;
                        kl(ii,jj)=trace(pinv(cj)*ci)-log(det(ci)/det(cj))-size(cj,1);
                    end;
                end;
            end; 

            first_level_results.pw_oasm(indind)=sum(squash(kl))/(NK^2-NK);

            %%
        end;

        %%%%%%
        %% search for dir
        w=weights{1}';   

        if(size(w,2)>1)
        else
            oat.first_level.connectivity_seed_search_direction=0;
        end;

        tmp=w'*w;   

        [u, ~] = svd(tmp,'econ');                              
        yeta1 = (u(:,1));

        if(~oat.first_level.connectivity_seed_search_direction),
            yeta2 = zeros(size(yeta1));
        else
            yeta2 = (u(:,2));
        end;

        w=weights_seed{1}';                
        tmp=w'*w;   

        [u, ~] = svd(tmp,'econ');                                  
        xeta1 = (u(:,1));

        if(~oat.first_level.connectivity_seed_search_direction),
            xeta2 = zeros(size(xeta1));
        else
            xeta2 = (u(:,2));
        end;

        if(oat.first_level.connectivity_seed_search_direction),
            as=0:0.4:pi;                
            bs=0:0.4:pi;               
        else
            as=0;
            bs=0;
        end;
        
        if(indind>0),

            plv=zeros(length(as),length(bs));
            he=zeros(length(as),length(bs));

            for ff=1:size(conn.dat_seed_tf,2),
                
            clear cope varcope;
            for a=1:length(as),

                % now do windowing over time                    
                x=permute(conn.dat_seed_tf(:,ff,1,:),[1 4 2 3]);
                xin=(x*(xeta1*cos(as(a))+xeta2*sin(as(a))));

                % now compute HE:
                x=(sqrt(xin.*conj(xin)));

                % moving average 
                if ~isempty(first_level.time_moving_av_win_size)

                    windowSize = first_level.time_moving_av_win_size*fsample_posttf;

                    x=moving(x,round(windowSize));

                end
                
                % temporally downsample
                if post_movingaverage_ds_factor~=1,
                    x = osl_spm_resample(x',post_movingaverage_ds_factor)';  
                    
                    if(indind==1 && ff==1),
                        disp(['Sample rate data post moving average is ' num2str(fsample_posttf*post_movingaverage_ds_factor) ' Hz']);
                    end;
                end;
                        
                x=normalise(x);
                pinvxtx=pinv(x'*x);
                pinvx=pinv(x);

                x_phase=angle(xin);

                xraw=permute(conn.dat_seed,[1 4 2 3]);
                xraw=(xraw*(xeta1*cos(as(a))+xeta2*sin(as(a))));

                for b=1:length(bs),

                    if any(strcmp(oat.first_level.connectivity_fc_methods,'he')) || any(strcmp(oat.first_level.connectivity_fc_methods,'plv')),
                        
                        % now do windowing over time                    
                        y=permute(conn.dat_tf(:,ff,1,:),[1 4 2 3]);
                
                        yin=(y*(yeta1*cos(bs(b))+yeta2*sin(bs(b))));

                        % now compute HE:
                        y=(sqrt(yin.*conj(yin)));

                        % moving average 
                        if ~isempty(first_level.time_moving_av_win_size)

                            windowSize = first_level.time_moving_av_win_size*fsample_posttf;

                            y=moving(y,round(windowSize));

                        end

                        % temporally downsample
                        if post_movingaverage_ds_factor~=1,
                            y = osl_spm_resample(y',post_movingaverage_ds_factor)';                  
                        end;
               
                        y=normalise(y);

                        y_phase=angle(yin);

                        [cope, varcope]=glm_fast_for_meg(y,x,pinvxtx,pinvx,[1]);

                        he(a,b)=cope/sqrt(varcope);

                        plv(a,b) = fisher_xform(abs(mean(exp(sqrt(-1)*(y_phase-x_phase))))); % Measuring Phase Synchrony in Brain Signals. Jean-Philippe Lachaux et al, HBM 1999
                    end;

                    if any(strcmp(oat.first_level.connectivity_fc_methods,'synchro')),
                        M=10;
                        deltasecs=1/(M*bandcentre);
                        delta=ceil(deltasecs/tres);

                        yraw=permute(conn.dat,[1 4 2 3]);
                
                        yraw=(yraw*(yeta1*cos(bs(b))+yeta2*sin(bs(b))));
                        gs(a,b)=synchro([normalise(yraw) normalise(xraw)],M,delta,10,10);
                    end;
                end;

            end;

            tmp=he;

            [i,j]=find(tmp==max(max(tmp)));                
            first_level_results.he_angles_seed(indind,ff)=as(i(1));                
            first_level_results.hes(indind,ff)=max(max(tmp));
        
            first_level_results.cope(indind,ff)=first_level_results.hes(indind,ff);       
            first_level_results.stdcope(indind,ff)=1;       

            tmp=plv;

            
            [i,j]=find(tmp==max(max(tmp)));    

            first_level_results.plv_angles_seed(indind,ff)=as(i(1));                
            first_level_results.plvs(indind,ff)=max(max(tmp));
            
        end; %end ff freq for loop
        
            if any(strcmp(oat.first_level.connectivity_fc_methods,'synchro')),

                tmp=gs;

                [i,j]=find(tmp==max(max(tmp)));                
                first_level_results.gs_angles_seed(indind)=as(i);                
                first_level_results.gss(indind)=max(max(tmp));
            
            end;
            
        end;           
         
    end; % for indind = 1:Nvoxels_out
    end; % if any(strcmp(oat.first_level.connectivity_fc_methods
    
    if any(strcmp(oat.first_level.connectivity_fc_methods,'hmm_kl')),
    
        disp('Doing HMM...');           
                 
        NK=6;
        hmm_num_starts=1;
        hmm_kl_use_source_seed=0;
        hmm_pca_dim=40;
        use_seed=1;
        
        disp(['Inferring ' num2str(NK) ' HMM states']);            
        
        Sh=[];                

        if(hmm_kl_use_source_seed)
            sensor_data_zerolagcorr_tmp=permute(sensor_data_zerolagcorr,[2 1 3]);
            %sensor_data_zerolagcorr_tmp=permute(sensor_data,[2 1 3]);
            sd=[size(sensor_data_zerolagcorr_tmp) 1 1]; 
            sensor_data_zerolagcorr_tmp=reshape(sensor_data_zerolagcorr_tmp,[sd(1), sd(2)*sd(3)]);
            Sh.data=sensor_data_zerolagcorr_tmp;
            Sh.data=[Sh.data, conn.dat_seed];
            
            if(hmm_pca_dim>0)
                disp('Doing pre-HMM PCA...');
                [allsvd,Apca]=pca(Sh.data,hmm_pca_dim);
                pinvApca=pinv(Apca);              
            else
                error('hmm_pca_dim unspecified');
            end;
            
            Sh.data=(pinvApca*Sh.data')';    
        else,
            sensor_data_tmp=permute(sensor_data,[2 1 3]);
            %sensor_data_zerolagcorr_tmp=permute(sensor_data,[2 1 3]);
            sd=[size(sensor_data_tmp) 1 1]; 
            sensor_data_tmp=reshape(sensor_data_tmp,[sd(1), sd(2)*sd(3)]);
            Sh.data=sensor_data_tmp;
            
            if(hmm_pca_dim>0)
                disp('Doing pre-HMM PCA...');
                [allsvd,Apca]=pca(Sh.data,hmm_pca_dim);
                pinvApca=pinv(Apca);              
            else
                error('hmm_pca_dim unspecified');
            end;
            
            Sh.data=(pinvApca*Sh.data')';    
        end;
        
        % add seed onto the end:
        
        Sh.data=normalise(Sh.data);
        
        Sh.NK=NK;

        Sh.num_starts=hmm_num_starts;

        [ hmm, block ]=run_multistart_hmm(Sh);                
        
        NK=hmm.K;       
        
        Shh=[];
        Shh.NK=NK;
        Shh.hmm=hmm;        
        Shh.block=block;
        Shh.tres=1/fsample_tf;  
        Shh.Apca=Apca;
        
        plot_hmm( Shh );
  
        cv=Apca*hmm.state(1).Cov*Apca';
        
        if(hmm_kl_use_source_seed)
            weights_seed=zeros(2,size(cv,1));
            weights_seed(:,size(cv,1)-1:end)=eye(2);    
        else
            weights_seed=source_recon_results.BF.inverse.W.MEG{1}{first_level_results.seed_index}; % needs to be changed to get working with HMM beamformer
            tmp=weights_seed*weights_seed';   
            [u, ~] = svd(tmp,'econ');                              
            weights_seed = u(:,1:2)'*weights_seed;
            weights_seed=weights_seed/sqrt(trace(weights_seed*weights_seed'));
            weights_seed=[weights_seed];
        end;
                
        if(numel(source_recon_results.BF.inverse.W.MEG)>1),
            error('Not implemented for HMM beamformer');
        end;
        
        if(use_seed)
            kl=zeros(Nvoxels_out,NK,NK);
        
            ft_progress('init', 'etf');
            for voxind=1:Nvoxels_out;

                ft_progress(voxind/Nvoxels_out);                                   

                clear cv_pair;
                for kk=1:NK,

                    cv=Apca*hmm.state(kk).Cov*Apca';

                    weights=source_recon_results.BF.inverse.W.MEG{1}{voxind}; % needs to be changed to get working with HMM beamformer
                    tmp=weights*weights';   
                    [u, ~] = svd(tmp,'econ');                              
                    weights = u(:,1:2)'*weights;
                    weights=weights/sqrt(trace(weights*weights'));

                    if(hmm_kl_use_source_seed)
                        weights=[weights, zeros(2,2)];
                    end;

                    weights_pair=[weights', weights_seed']';

                    cv_pair{kk}=weights_pair*cv*weights_pair';   

                end;            

               for jj=1:NK,
                    cj=cv_pair{jj};
                    pcj=inv(cj);
                    ldcj_mNj=logdet(cj,'chol')-size(cj,1);

                    for ii=jj+1:NK,   
                    %for jj=1:NK,
                    %    if(ii~=jj),    
                        ci=cv_pair{ii};

                        %  kl(voxind,ii,jj)=trace(pinv(cj)*ci)-log(det(ci)/det(cj))-size(cj,1);
                        kl(voxind,ii,jj)=trace(pcj*ci)-logdet(ci,'chol')+ldcj_mNj;  
                    %end;
                    end;
                end;   
                                                                        
            end;  
            
            ft_progress('close');
        else,
                    
            % precompute weights
            for voxind=1:Nvoxels_out;
                
                weights=source_recon_results.BF.inverse.W.MEG{1}{voxind}; % needs to be changed to get working with HMM beamformer
                tmp=weights*weights';   
                [u, ~] = svd(tmp,'econ');                              
                weights = u(:,1:2)'*weights;
                weights=weights/sqrt(trace(weights*weights'));
                weights_stored{voxind}=weights;
            end;
            
            % precompute state cvs:
            for kk=1:NK,
                cvs_stored{kk}=Apca*hmm.state(kk).Cov*Apca';
            end;
            
            kl=zeros(Nvoxels_out,Nvoxels_out,NK,NK);
        
            ft_progress('init', 'etf');
            
            cv_pair=zeros(size(weights_stored{1},1)*2,size(weights_stored{1},1)*2,NK);

            for voxind_seed=1:Nvoxels_out;

                ft_progress(voxind_seed/Nvoxels_out);  

                for voxind=1:Nvoxels_out;

                    if(voxind_seed~=voxind),                                                                     

                        clear cv_pair;
                        for kk=1:NK,

                            weights_pair=[weights_stored{voxind}', weights_stored{voxind_seed}']';
                            cv_pair(:,:,kk)=weights_pair*cvs_stored{kk}*weights_pair';   

                        end;            


                        for jj=1:NK,
                            cj=cv_pair(:,:,jj);
                            pcj=inv(cj);
                            ldcj_mNj=logdet(cj,'chol')-size(cj,1);

                            for ii=jj+1:NK,     
                                ci=cv_pair(:,:,ii);

                                %  kl(voxind,ii,jj)=trace(pinv(cj)*ci)-log(det(ci)/det(cj))-size(cj,1);
                                kl(voxind_seed,voxind,ii,jj)=trace(pcj*ci)-logdet(ci,'chol')+ldcj_mNj;  
                                
                            end;
                        end;        
                    end;
                end;
            end;
           
            ft_progress('close');
        end;
                
        %% node OASM
        % cd ~; load kl;
        
        first_level_results.kl_edge_oasm=sum(sum(kl,3),4);
        
        first_level_results.kl_node_oasm=sum(sum(sum(kl,2),3),4);
   
        first_level_results.klm_edge_oasm=max(max(kl,[],3),[],4);
        
        first_level_results.klm_node_oasm=max(max(max(kl,[],2),[],3),[],4);
      
        first_level_results.kl_edge_index=(1:size(first_level_results.kl_edge_oasm,1))';
    end;
   
                
    first_level_results.source_recon_results_fname=source_recon_results_fname;
    first_level_results.first_level=first_level;
    first_level_results.source_recon=source_recon_results.source_recon;
    first_level_results.level=1;
    %first_level_results.source_recon_results_fname=source_recon_results.fname;
    first_level_results.name=first_level.name;
    first_level_results.recon_method=source_recon_results.recon_method;
    first_level_results.session_name=source_recon_results.session_name;
          
    first_level_results.fname=[ first_level_results.session_name '_' first_level_results.name ];
    disp(['Saving statistics in file ' oat.source_recon.dirname '/' first_level_results.fname]);
    osl_save_oat_results( oat, first_level_results );
    results_fnames{sub}=first_level_results.fname;

    if (strcmp(source_recon_results.recon_method,'none')),
    else
        disp(['To create niftii files from this use a call to osl_save_nii_stats']);
    end;
    
end;

function y=fisher_xform(x)

y=0.5*log((1+x)/(1-x));
    

